Molecular signatures underlying heterogeneous hypertrophic responsiveness to resistance trainin in older men and women: a within-subject design

Author

Manoel Lixandrao

Published

August 2, 2025

Install and load libraries

# Packages from the Bioconductor
if(!require(BiocManager)) install.packages("BiocManager")
if(!require(EnhancedVolcano)) BiocManager::install("EnhancedVolcano")
if(!require(limma)) BiocManager::install("limma")
if(!require(Biobase)) BiocManager::install("Biobase")
if(!require(edgeR)) BiocManager::install("edgeR")
if(!require(biomaRt)) BiocManager::install("biomaRt")
if(!require(org.Hs.eg.db)) BiocManager::install("org.Hs.eg.db")
if(!require(DESeq2)) BiocManager::install("DESeq2")
if(!require(topGO)) BiocManager::install("topGO")
if(!require(RISmed)) BiocManager::install("RISmed")
if(!require(pathview)) BiocManager::install("pathview")
if(!require(gage)) BiocManager::install("gage")
if(!require(gageData)) BiocManager::install("gageData")
if(!require(ReactomePA)) BiocManager::install("ReactomePA")
if(!require(clusterProfiler)) BiocManager::install("clusterProfiler")
if(!require(enrichplot)) BiocManager::install("enrichplot")
if(!require(DOSE)) BiocManager::install("DOSE")

# Packages from the github
if(!require(PLIER)) devtools::install_github("wgmao/PLIER")

# Packages from the r
if(!require(ggpubr)) install.packages("ggpubr")
if(!require(ggsignif)) install.packages("ggsignif")
if(!require(matrixStats)) install.packages("matrixStats")
if(!require(matrixStats)) install.packages("matrixStats")
if(!require(RColorBrewer)) install.packages("RColorBrewer")
if(!require(gplots)) install.packages("gplots")
if(!require(ggplot2)) install.packages("ggplot2")
if(!require(dplyr)) install.packages("dplyr")
if(!require(tidyverse)) install.packages("tidyverse")
if(!require(readxl)) install.packages("readxl")
if(!require(rstatix)) install.packages("rstatix")
if(!require(kableExtra)) install.packages("kableExtra")
if(!require(broom)) install.packages("broom")
if(!require(emmeans)) install.packages("emmeans")
if(!require(lme4)) install.packages("lme4")
if(!require(easyPubMed)) install.packages("easyPubMed")
if(!require(factoextra)) install.packages("factoextra")
if(!require(gt)) install.packages("gt")
if(!require(cowplot)) install.packages("cowplot")

Data

Load data

# List files and load into R environment
list.files("0_data", full.names = TRUE, pattern = ".Rds") |>
  walk(\(x) load(x, envir = .GlobalEnv))

# Print objects
gt::gt(design.data)
SubjectID Sex Group TimePoint Set 1_set 4_set condition dummyId
A_102 1 L1_L4 pre 1 D E baseline 1
A_104 1 H1_H4 pre 1 E D baseline 2
A_105 2 L1_H4 pre 1 E D baseline 3
A_106 2 L1_H4 pre 1 E D baseline 4
A_107 1 L1_H4 pre 1 E D baseline 5
A_110 2 L1_H4 pre 1 D E baseline 6
A_114 2 L1_H4 pre 1 E D baseline 7
A_12 1 L1_L4 pre 1 E D baseline 8
A_17 1 L1_H4 pre 1 D E baseline 9
A_27 1 H1_H4 pre 1 E D baseline 10
A_30 2 L1_L4 pre 1 D E baseline 11
A_32 1 H1_H4 pre 1 E D baseline 12
A_35 2 L1_L4 pre 1 E D baseline 13
A_36 2 H1_H4 pre 1 D E baseline 14
A_38 2 H1_H4 pre 1 E D baseline 15
A_50 2 L1_H4 pre 1 D E baseline 16
A_59 2 L1_L4 pre 1 D E baseline 17
A_5 1 L1_L4 pre 1 E D baseline 18
A_6 2 H1_H4 pre 1 E D baseline 19
A_76 1 L1_H4 pre 1 E D baseline 20
A_77 1 L1_L4 pre 1 E D baseline 21
A_80 1 H1_H4 pre 1 D E baseline 22
A_86 2 H1_H4 pre 1 E D baseline 23
A_90 1 L1_L4 pre 1 E D baseline 24
A_94 2 L1_L4 pre 1 D E baseline 25
A_97 2 H1_H4 pre 1 D E baseline 26
A_98 1 L1_H4 pre 1 E D baseline 27
BD_102 1 L1_L4 post 1 D E post_1set 1
BD_104 1 H1_H4 post 4 E D post_4set 2
BD_105 2 L1_H4 post 4 E D post_4set 3
BD_106 2 L1_H4 post 4 E D post_4set 4
BD_107 1 L1_H4 post 4 E D post_4set 5
BD_110 2 L1_H4 post 1 D E post_1set 6
BD_114 2 L1_H4 post 4 E D post_4set 7
BD_12 1 L1_L4 post 4 E D post_4set 8
BD_17 1 L1_H4 post 1 D E post_1set 9
BD_27 1 H1_H4 post 4 E D post_4set 10
BD_30 2 L1_L4 post 1 D E post_1set 11
BD_32 1 H1_H4 post 4 E D post_4set 12
BD_35 2 L1_L4 post 4 E D post_4set 13
BD_36 2 H1_H4 post 1 D E post_1set 14
BD_38 2 H1_H4 post 4 E D post_4set 15
BD_50 2 L1_H4 post 1 D E post_1set 16
BD_59 2 L1_L4 post 1 D E post_1set 17
BD_5 1 L1_L4 post 4 E D post_4set 18
BD_6 2 H1_H4 post 4 E D post_4set 19
BD_76 1 L1_H4 post 4 E D post_4set 20
BD_77 1 L1_L4 post 4 E D post_4set 21
BD_80 1 H1_H4 post 1 D E post_1set 22
BD_86 2 H1_H4 post 4 E D post_4set 23
BD_90 1 L1_L4 post 4 E D post_4set 24
BD_94 2 L1_L4 post 1 D E post_1set 25
BD_97 2 H1_H4 post 1 D E post_1set 26
BD_98 1 L1_H4 post 4 E D post_4set 27
BE_102 1 L1_L4 post 4 D E post_4set 1
BE_104 1 H1_H4 post 1 E D post_1set 2
BE_105 2 L1_H4 post 1 E D post_1set 3
BE_106 2 L1_H4 post 1 E D post_1set 4
BE_107 1 L1_H4 post 1 E D post_1set 5
BE_110 2 L1_H4 post 4 D E post_4set 6
BE_114 2 L1_H4 post 1 E D post_1set 7
BE_12 1 L1_L4 post 1 E D post_1set 8
BE_17 1 L1_H4 post 4 D E post_4set 9
BE_27 1 H1_H4 post 1 E D post_1set 10
BE_30 2 L1_L4 post 4 D E post_4set 11
BE_32 1 H1_H4 post 1 E D post_1set 12
BE_35 2 L1_L4 post 1 E D post_1set 13
BE_36 2 H1_H4 post 4 D E post_4set 14
BE_38 2 H1_H4 post 1 E D post_1set 15
BE_50 2 L1_H4 post 4 D E post_4set 16
BE_59 2 L1_L4 post 4 D E post_4set 17
BE_5 1 L1_L4 post 1 E D post_1set 18
BE_6 2 H1_H4 post 1 E D post_1set 19
BE_76 1 L1_H4 post 1 E D post_1set 20
BE_77 1 L1_L4 post 1 E D post_1set 21
BE_80 1 H1_H4 post 4 D E post_4set 22
BE_86 2 H1_H4 post 1 E D post_1set 23
BE_90 1 L1_L4 post 1 E D post_1set 24
BE_94 2 L1_L4 post 4 D E post_4set 25
BE_97 2 H1_H4 post 4 D E post_4set 26
BE_98 1 L1_H4 post 1 E D post_1set 27 
gt::gt(head(gene.list, 20))
ensembl_gene_id gene_biotype hgnc_symbol chromosome_name start_position end_position external_gene_name
ENSG00000210049 Mt_tRNA MT-TF MT 577 647 MT-TF
ENSG00000211459 Mt_rRNA MT-RNR1 MT 648 1601 MT-RNR1
ENSG00000210077 Mt_tRNA MT-TV MT 1602 1670 MT-TV
ENSG00000210082 Mt_rRNA MT-RNR2 MT 1671 3229 MT-RNR2
ENSG00000209082 Mt_tRNA MT-TL1 MT 3230 3304 MT-TL1
ENSG00000198888 protein_coding MT-ND1 MT 3307 4262 MT-ND1
ENSG00000210100 Mt_tRNA MT-TI MT 4263 4331 MT-TI
ENSG00000210107 Mt_tRNA MT-TQ MT 4329 4400 MT-TQ
ENSG00000210112 Mt_tRNA MT-TM MT 4402 4469 MT-TM
ENSG00000198763 protein_coding MT-ND2 MT 4470 5511 MT-ND2
ENSG00000210117 Mt_tRNA MT-TW MT 5512 5579 MT-TW
ENSG00000210127 Mt_tRNA MT-TA MT 5587 5655 MT-TA
ENSG00000210135 Mt_tRNA MT-TN MT 5657 5729 MT-TN
ENSG00000210140 Mt_tRNA MT-TC MT 5761 5826 MT-TC
ENSG00000210144 Mt_tRNA MT-TY MT 5826 5891 MT-TY
ENSG00000198804 protein_coding MT-CO1 MT 5904 7445 MT-CO1
ENSG00000210151 Mt_tRNA MT-TS1 MT 7446 7514 MT-TS1
ENSG00000210154 Mt_tRNA MT-TD MT 7518 7585 MT-TD
ENSG00000198712 protein_coding MT-CO2 MT 7586 8269 MT-CO2
ENSG00000210156 Mt_tRNA MT-TK MT 8295 8364 MT-TK 
expression.data |> 
  rownames_to_column("ENSEMBL") |> 
  (\(x) gt::gt(head(x, 20)))()
ENSEMBL A_102 A_104 A_105 A_106 A_107 A_110 A_114 A_12 A_17 A_27 A_30 A_32 A_35 A_36 A_38 A_50 A_59 A_5 A_6 A_76 A_77 A_80 A_86 A_90 A_94 A_97 A_98 BD_102 BD_104 BD_105 BD_106 BD_107 BD_110 BD_114 BD_12 BD_17 BD_27 BD_30 BD_32 BD_35 BD_36 BD_38 BD_50 BD_59 BD_5 BD_6 BD_76 BD_77 BD_80 BD_86 BD_90 BD_94 BD_97 BD_98 BE_102 BE_104 BE_105 BE_106 BE_107 BE_110 BE_114 BE_12 BE_17 BE_27 BE_30 BE_32 BE_35 BE_36 BE_38 BE_50 BE_59 BE_5 BE_6 BE_76 BE_77 BE_80 BE_86 BE_90 BE_94 BE_97 BE_98
ENSG00000223972 5 0 4 10 7 5 0 3 9 3 7 3 8 7 2 6 4 9 0 8 13 4 2 11 5 5 11 10 0 8 4 8 8 2 2 15 7 0 5 2 6 1 1 5 18 7 10 22 3 3 3 5 11 12 7 1 3 10 14 12 0 0 15 0 10 7 5 4 0 3 7 4 0 8 20 9 3 11 6 3 13
ENSG00000227232 612 811 746 618 695 857 613 558 910 605 854 1000 833 538 629 901 853 1532 732 1109 1064 479 522 844 580 616 517 562 739 679 559 945 728 639 538 982 674 774 905 680 679 469 982 1223 1418 962 1126 1044 516 448 801 442 440 554 527 761 536 709 857 846 769 608 1209 804 818 908 693 680 606 739 940 1249 774 1001 1285 433 451 610 594 628 628
ENSG00000278267 5 8 3 4 8 16 5 10 10 13 11 15 4 5 15 19 23 19 6 12 16 7 14 20 7 8 1 3 5 3 5 25 0 12 11 11 1 1 8 9 8 8 16 11 18 9 18 14 4 15 24 4 7 13 7 15 4 2 7 15 22 6 14 12 19 23 20 16 4 13 12 17 16 10 15 0 6 6 9 13 7
ENSG00000243485 1 10 14 1 4 2 1 6 3 5 23 12 17 3 5 0 10 2 0 4 0 8 4 3 4 4 4 2 3 6 0 1 24 1 6 3 4 23 5 4 0 4 3 9 11 5 1 4 8 4 2 4 5 1 1 2 1 8 2 11 0 3 10 13 7 14 14 13 8 2 7 1 5 0 3 4 2 0 6 9 1
ENSG00000274890 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
ENSG00000237613 2 0 0 0 0 1 0 0 1 0 1 1 0 0 0 0 0 0 0 4 0 0 0 1 0 0 1 0 1 2 0 0 0 0 0 0 1 0 3 2 0 2 0 1 1 1 1 1 1 0 2 0 0 3 0 0 0 1 0 0 0 0 0 0 1 5 0 0 1 0 2 0 0 3 2 1 0 0 0 1 0
ENSG00000268020 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
ENSG00000240361 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
ENSG00000186092 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
ENSG00000238009 11 3 13 7 2 15 18 3 18 1 10 16 15 2 0 11 19 4 6 10 7 15 5 13 4 4 17 25 1 17 12 5 17 2 26 9 15 5 5 5 0 2 14 18 4 1 12 18 11 3 21 11 3 14 6 5 13 7 5 14 12 7 9 6 9 3 7 3 0 15 5 6 1 13 5 8 14 23 11 7 12
ENSG00000239945 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
ENSG00000233750 10 7 5 1 16 14 9 3 7 3 23 10 12 5 7 15 9 8 5 11 7 6 6 11 12 10 4 11 5 9 9 14 17 5 5 10 3 11 11 7 0 6 10 9 6 11 4 15 13 6 9 10 5 3 14 8 5 16 4 10 13 22 13 24 16 9 12 7 13 24 8 4 4 7 9 14 9 5 6 3 12
ENSG00000268903 43 46 51 39 39 104 38 58 51 41 67 75 17 34 30 17 50 90 50 50 41 32 29 33 55 30 32 50 40 41 33 49 128 18 55 56 40 43 73 21 50 23 34 61 75 57 22 45 23 32 54 28 33 33 63 41 28 45 60 78 36 31 79 42 79 76 24 57 49 25 33 63 44 29 53 43 33 34 36 40 59
ENSG00000269981 18 14 16 13 32 34 4 17 22 14 16 8 8 24 13 5 15 44 35 7 8 15 14 11 13 14 11 25 8 8 8 34 62 10 31 20 9 16 11 4 35 18 12 35 55 56 14 38 7 7 6 9 3 3 13 24 13 6 13 23 13 21 19 30 15 11 11 33 7 5 13 18 23 11 30 7 21 10 28 16 9
ENSG00000239906 6 9 11 7 11 14 1 7 11 18 14 7 2 12 6 5 5 8 8 6 4 4 8 4 4 4 0 9 9 8 5 10 5 5 6 12 8 4 2 6 10 16 2 3 27 8 3 2 2 30 5 2 3 3 19 13 5 6 7 8 6 4 13 3 2 8 1 15 11 1 4 11 8 15 8 2 31 5 13 4 5
ENSG00000241860 201 95 246 115 178 226 131 153 160 158 154 133 103 183 131 125 153 203 157 136 136 139 119 132 118 112 112 204 77 211 97 334 171 97 137 192 163 174 121 88 304 138 115 164 238 256 135 198 121 118 86 112 60 126 211 98 187 93 285 217 108 139 223 157 130 98 103 320 151 106 128 193 241 142 163 94 108 100 134 92 146
ENSG00000222623 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 2 0 1 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
ENSG00000241599 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0
ENSG00000279928 2 5 1 3 14 5 3 0 9 0 6 4 11 5 0 4 1 3 0 8 10 1 0 9 7 2 8 6 0 11 2 7 13 0 1 4 2 2 5 2 6 1 0 4 4 5 5 17 2 1 3 3 7 6 11 2 0 3 17 8 3 1 4 0 11 4 5 1 0 0 0 4 0 3 15 2 2 7 3 1 10
ENSG00000279457 639 698 736 604 601 951 887 663 882 763 811 978 805 565 716 1039 999 1485 802 955 1060 607 489 933 704 802 618 557 632 625 523 917 908 695 687 1009 680 666 783 701 623 557 1165 1370 1332 1038 899 860 563 506 910 536 525 612 578 640 479 690 817 985 1113 632 1186 980 759 866 765 757 553 967 1139 1190 867 879 1383 516 579 685 592 723 664

Data wrangling

# Merge TimePoint and Set into Condition (baseline, post_1set, post_4set) for the RNA sequencing analysis.
design.data <- design.data |> 
  mutate(condition = paste(TimePoint, Set, sep = "_"),
         condition = factor(condition, levels = c("pre_1", "post_1", "post_4"), labels = c("baseline", "post_1set", "post_4set"))) |> 
  # create a dummyID to input in the voom-limma pipeline (within subject correlation)
  mutate(dummyId = rep(1:27, 3))

# Set rownames for your design.data 
rownames(design.data) <- design.data$SubjectID

gt::gt(design.data)
SubjectID Sex Group TimePoint Set 1_set 4_set condition dummyId
A_102 1 L1_L4 pre 1 D E baseline 1
A_104 1 H1_H4 pre 1 E D baseline 2
A_105 2 L1_H4 pre 1 E D baseline 3
A_106 2 L1_H4 pre 1 E D baseline 4
A_107 1 L1_H4 pre 1 E D baseline 5
A_110 2 L1_H4 pre 1 D E baseline 6
A_114 2 L1_H4 pre 1 E D baseline 7
A_12 1 L1_L4 pre 1 E D baseline 8
A_17 1 L1_H4 pre 1 D E baseline 9
A_27 1 H1_H4 pre 1 E D baseline 10
A_30 2 L1_L4 pre 1 D E baseline 11
A_32 1 H1_H4 pre 1 E D baseline 12
A_35 2 L1_L4 pre 1 E D baseline 13
A_36 2 H1_H4 pre 1 D E baseline 14
A_38 2 H1_H4 pre 1 E D baseline 15
A_50 2 L1_H4 pre 1 D E baseline 16
A_59 2 L1_L4 pre 1 D E baseline 17
A_5 1 L1_L4 pre 1 E D baseline 18
A_6 2 H1_H4 pre 1 E D baseline 19
A_76 1 L1_H4 pre 1 E D baseline 20
A_77 1 L1_L4 pre 1 E D baseline 21
A_80 1 H1_H4 pre 1 D E baseline 22
A_86 2 H1_H4 pre 1 E D baseline 23
A_90 1 L1_L4 pre 1 E D baseline 24
A_94 2 L1_L4 pre 1 D E baseline 25
A_97 2 H1_H4 pre 1 D E baseline 26
A_98 1 L1_H4 pre 1 E D baseline 27
BD_102 1 L1_L4 post 1 D E post_1set 1
BD_104 1 H1_H4 post 4 E D post_4set 2
BD_105 2 L1_H4 post 4 E D post_4set 3
BD_106 2 L1_H4 post 4 E D post_4set 4
BD_107 1 L1_H4 post 4 E D post_4set 5
BD_110 2 L1_H4 post 1 D E post_1set 6
BD_114 2 L1_H4 post 4 E D post_4set 7
BD_12 1 L1_L4 post 4 E D post_4set 8
BD_17 1 L1_H4 post 1 D E post_1set 9
BD_27 1 H1_H4 post 4 E D post_4set 10
BD_30 2 L1_L4 post 1 D E post_1set 11
BD_32 1 H1_H4 post 4 E D post_4set 12
BD_35 2 L1_L4 post 4 E D post_4set 13
BD_36 2 H1_H4 post 1 D E post_1set 14
BD_38 2 H1_H4 post 4 E D post_4set 15
BD_50 2 L1_H4 post 1 D E post_1set 16
BD_59 2 L1_L4 post 1 D E post_1set 17
BD_5 1 L1_L4 post 4 E D post_4set 18
BD_6 2 H1_H4 post 4 E D post_4set 19
BD_76 1 L1_H4 post 4 E D post_4set 20
BD_77 1 L1_L4 post 4 E D post_4set 21
BD_80 1 H1_H4 post 1 D E post_1set 22
BD_86 2 H1_H4 post 4 E D post_4set 23
BD_90 1 L1_L4 post 4 E D post_4set 24
BD_94 2 L1_L4 post 1 D E post_1set 25
BD_97 2 H1_H4 post 1 D E post_1set 26
BD_98 1 L1_H4 post 4 E D post_4set 27
BE_102 1 L1_L4 post 4 D E post_4set 1
BE_104 1 H1_H4 post 1 E D post_1set 2
BE_105 2 L1_H4 post 1 E D post_1set 3
BE_106 2 L1_H4 post 1 E D post_1set 4
BE_107 1 L1_H4 post 1 E D post_1set 5
BE_110 2 L1_H4 post 4 D E post_4set 6
BE_114 2 L1_H4 post 1 E D post_1set 7
BE_12 1 L1_L4 post 1 E D post_1set 8
BE_17 1 L1_H4 post 4 D E post_4set 9
BE_27 1 H1_H4 post 1 E D post_1set 10
BE_30 2 L1_L4 post 4 D E post_4set 11
BE_32 1 H1_H4 post 1 E D post_1set 12
BE_35 2 L1_L4 post 1 E D post_1set 13
BE_36 2 H1_H4 post 4 D E post_4set 14
BE_38 2 H1_H4 post 1 E D post_1set 15
BE_50 2 L1_H4 post 4 D E post_4set 16
BE_59 2 L1_L4 post 4 D E post_4set 17
BE_5 1 L1_L4 post 1 E D post_1set 18
BE_6 2 H1_H4 post 1 E D post_1set 19
BE_76 1 L1_H4 post 1 E D post_1set 20
BE_77 1 L1_L4 post 1 E D post_1set 21
BE_80 1 H1_H4 post 4 D E post_4set 22
BE_86 2 H1_H4 post 1 E D post_1set 23
BE_90 1 L1_L4 post 1 E D post_1set 24
BE_94 2 L1_L4 post 4 D E post_4set 25
BE_97 2 H1_H4 post 4 D E post_4set 26
BE_98 1 L1_H4 post 1 E D post_1set 27

Data check

# Test if all subjects are in the expression.data. Should return TRUE
all(rownames(design.data) == colnames(expression.data))
[1] TRUE

RNA data (expression matrix)

Read in count-data and organize sample information

# Read count-data
x <- DGEList(counts = expression.data)

# Organize sample information
x$samples$sex <- factor(design.data$Sex[match(rownames(x$samples), rownames(design.data))])
x$samples$group <- factor(design.data$Group[match(rownames(x$samples), rownames(design.data))])
x$samples$timepoint <- factor(design.data$TimePoint[match(rownames(x$samples), rownames(design.data))])
x$samples$condition <- factor(design.data$condition[match(rownames(x$samples), rownames(design.data))])
x$samples$dummyid <- factor(design.data$dummyId[match(rownames(x$samples), rownames(design.data))])

# Organize gene list and names
x$genes <- gene.list[match(rownames(x), rownames(gene.list)),]

Remove low-expressed genes

keep.exprs <- filterByExpr(x, group = x$sample$group)
x <- x[keep.exprs, , keep.lib.sizes=FALSE]

# Check the size of your dataset 
dim(x)
[1] 16617    81

Normalizing gene expression distributions

x <- calcNormFactors(x, method = "TMM")

Annotate to the Ensembl

# Create a vector with correspond ENSEMBL to ENTREZID - add to the x$genes
x$genes["ENTREZID"] = mapIds(org.Hs.eg.db,
                  keys = rownames(x),
                  column = "ENTREZID",
                  keytype = "ENSEMBL",
                  multiVals = "first")

gt::gt(head(x$genes, 20))
ensembl_gene_id gene_biotype hgnc_symbol chromosome_name start_position end_position external_gene_name ENTREZID
ENSG00000227232 unprocessed_pseudogene WASH7P 1 14404 29570 WASH7P NA
ENSG00000278267 miRNA MIR6859-1 1 17369 17436 MIR6859-1 102466751
ENSG00000238009 lncRNA 1 89295 133723 AL627309.1 NA
ENSG00000233750 processed_pseudogene CICP27 1 131025 134836 CICP27 100420257
ENSG00000268903 processed_pseudogene 1 135141 135895 AL627309.6 NA
ENSG00000269981 processed_pseudogene 1 137682 137965 AL627309.7 NA
ENSG00000241860 lncRNA 1 141474 173862 AL627309.5 NA
ENSG00000279457 unprocessed_pseudogene WASH9P 1 185217 195411 WASH9P NA
ENSG00000228463 transcribed_processed_pseudogene 1 257864 359681 AP006222.1 728481
ENSG00000237094 transcribed_unprocessed_pseudogene 1 365389 522928 AL732372.2 NA
ENSG00000250575 unprocessed_pseudogene 1 491225 493241 AL732372.3 NA
ENSG00000230021 transcribed_processed_pseudogene 1 586071 827796 AL669831.3 NA
ENSG00000225972 unprocessed_pseudogene MTND1P23 1 629062 629433 MTND1P23 100887749
ENSG00000225630 unprocessed_pseudogene MTND2P28 1 629640 630683 MTND2P28 100652939
NA NA NA NA NA NA NA NA
ENSG00000237973 unprocessed_pseudogene MTCO1P12 1 631074 632616 MTCO1P12 107075141
ENSG00000229344 unprocessed_pseudogene MTCO2P12 1 632757 633438 MTCO2P12 107075310
ENSG00000240409 unprocessed_pseudogene MTATP8P1 1 633535 633741 MTATP8P1 106480795
ENSG00000248527 unprocessed_pseudogene MTATP6P1 1 633696 634376 MTATP6P1 106480796
ENSG00000198744 unprocessed_pseudogene MTCO3P12 1 634376 634922 MTCO3P12 107075270

Voom-Limma

# Create a group variable combining factors (group and condition)
group <- factor(paste(x$samples$group, x$samples$condition, sep = "."))
sex <- factor(x$samples$sex)

# Set the model.matrix
design <- model.matrix(~ 0 + group + sex)
colnames(design)[1:9] <- levels(group)

# double voom to define the individuals (x$samples$dummyid) as random factor
## Initiate the voom
eset <- voom(x, design, plot = TRUE)

## calculate the within-subject correlation
corfit <- duplicateCorrelation(eset, design, block = x$samples$dummyid)

## run the voom again using the corfit correlation
eset <- voom(x, design, block = x$samples$dummyid, correlation = corfit$consensus.correlation, plot = TRUE)

## fit the linear model
fit <- lmFit(eset, design, block = x$samples$dummyid, correlation = corfit$consensus.correlation)

## Create a full mode contrast
full.mode.contrast <- makeContrasts(
  
  # within-group comparison low-responders
  llBasalVs1 = L1_L4.post_1set-L1_L4.baseline,
  llBasalVs4 = L1_L4.post_4set-L1_L4.baseline,
  ll1Vs4 = L1_L4.post_4set-L1_L4.post_1set,
  
  # within-group comparison medium-responders
  lhBasalVs1 = L1_H4.post_1set-L1_H4.baseline,
  lhBasalVs4 = L1_H4.post_4set-L1_H4.baseline,
  lh1Vs4 = L1_H4.post_4set-L1_H4.post_1set,
  
  # within-group comparison high-responders
  hhBasalVs1 = H1_H4.post_1set-H1_H4.baseline,
  hhBasalVs4 = H1_H4.post_4set-H1_H4.baseline,
  hh1Vs4 = H1_H4.post_4set-H1_H4.post_1set,
  
  # Between-group comparison
  ll1Vslh1 = L1_H4.post_1set-L1_L4.post_1set,
  ll1Vshh1 = H1_H4.post_1set-L1_L4.post_1set,
  lh1Vshh1 = H1_H4.post_1set-L1_H4.post_1set,
  ll4Vslh4 = L1_H4.post_4set-L1_L4.post_4set,
  ll4Vshh4 = H1_H4.post_4set-L1_L4.post_4set,
  lh4Vshh4 = H1_H4.post_4set-L1_H4.post_4set,
  levels = design)

## Initiate the contrast comparisons
tmp <- contrasts.fit(fit, full.mode.contrast)

## Run the comparisons
tmp <- eBayes(tmp)

## Print your contrast results
tmp_results <- decideTests(tmp, p.value = 0.05)
summary(tmp_results)
       llBasalVs1 llBasalVs4 ll1Vs4 lhBasalVs1 lhBasalVs4 lh1Vs4 hhBasalVs1
Down            0          0      0         60         22      0          0
NotSig      16617      16617  16617      16529      16585  16617      16617
Up              0          0      0         28         10      0          0
       hhBasalVs4 hh1Vs4 ll1Vslh1 ll1Vshh1 lh1Vshh1 ll4Vslh4 ll4Vshh4 lh4Vshh4
Down           22      0        0        0        2        0        0        0
NotSig      16577  16617    16617    16617    16614    16617    16617    16617
Up             18      0        0        0        1        0        0        0
## extract the log fold change for every comparison of interst
llBasalVs1 <- topTreat(tmp, coef = "llBasalVs1", n=Inf)
llBasalVs4 <- topTreat(tmp, coef = "llBasalVs4", n=Inf)
ll1Vs4 <- topTreat(tmp, coef = "ll1Vs4", n=Inf)

lhBasalVs1 <- topTreat(tmp, coef = "lhBasalVs1", n=Inf)
lhBasalVs4 <- topTreat(tmp, coef = "lhBasalVs4", n=Inf)
lh1Vs4 <- topTreat(tmp, coef = "lh1Vs4", n=Inf)

hhBasalVs1 <- topTreat(tmp, coef = "hhBasalVs1", n=Inf)
hhBasalVs4 <- topTreat(tmp, coef = "hhBasalVs4", n=Inf)
hh1Vs4 <- topTreat(tmp, coef = "hh1Vs4", n=Inf)

lh1Vshh1 <- topTreat(tmp, coef = "lh1Vshh1", n=Inf)

Filter DE genes

DE_list <- list(llBasalVs1, llBasalVs4, ll1Vs4, lhBasalVs1, lhBasalVs4, lh1Vs4, hhBasalVs1, hhBasalVs4, hh1Vs4) |> 
  map(\(x) as.data.frame(x)) |> 
  map(\(x) x |> mutate(gene_id = rownames(x)) |> 
      dplyr::select(-ensembl_gene_id) |> 
      dplyr::select("ENSEMBL ID" = gene_id, "Gene name" = external_gene_name, ENTREZID, "Log fold-change" = logFC, "P-value" = P.Value, "Adjusted P-value" = adj.P.Val)
      ) |> 
  set_names(list("llBasalVs1", "llBasalVs4", "ll1Vs4", "lhBasalVs1", "lhBasalVs4", "lh1Vs4", "hhBasalVs1", "hhBasalVs4", "hh1Vs4")) 

# Add labels (gene name) based on the p-values
DE_list_filtered <-  DE_list |> 
  map(\(x) x |> arrange(`P-value`) |> 
        mutate(labels = if_else(`Adjusted P-value` <= 0.05, `Gene name`, NA)))

groups_names <- paste(rep(c("Low-responders", "Medium-responders", "High-responders"), each = 3), c("Baseline vs 1 set", "Baseline vs 4 sets", "1 vs 4 sets"))

DE_list_filtered |>
  map(\(x) head(x, 20)) |>
  map2(groups_names, \(x, header) gt(x) |> tab_header(title = header))
$llBasalVs1
Low-responders Baseline vs 1 set
ENSEMBL ID Gene name ENTREZID Log fold-change P-value Adjusted P-value labels
ENSG00000204084 INPP5B 3633 -0.3623201 0.0001105966 0.8349042 NA
ENSG00000279434 AL049776.1 NA -0.9933954 0.0001316651 0.8349042 NA
ENSG00000101000 PROCR 10544 1.0043615 0.0001821190 0.8349042 NA
ENSG00000262728 AC123768.2 101927788 -1.3407008 0.0002917917 0.8349042 NA
ENSG00000267060 PTGES3L 100885848 0.3163270 0.0004875871 0.8349042 NA
ENSG00000176697 BDNF 627 -2.0823783 0.0005267215 0.8349042 NA
ENSG00000157985 AGAP1 116987 0.2705296 0.0006045699 0.8349042 NA
ENSG00000270728 AL035413.1 NA -0.5305436 0.0006383850 0.8349042 NA
ENSG00000196453 ZNF777 27153 0.2453884 0.0008027980 0.8349042 NA
ENSG00000170385 SLC30A1 7779 0.7387403 0.0008397693 0.8349042 NA
ENSG00000100150 DEPDC5 9681 0.3048450 0.0008586730 0.8349042 NA
ENSG00000131788 PIAS3 10401 -0.2375297 0.0009221041 0.8349042 NA
ENSG00000159200 RCAN1 1827 -0.6318909 0.0009518243 0.8349042 NA
ENSG00000143458 GABPB2 126626 -0.3029994 0.0009683269 0.8349042 NA
ENSG00000170153 RNF150 57484 0.5134830 0.0010010480 0.8349042 NA
ENSG00000169131 ZNF354A 6940 0.4165301 0.0011852225 0.8349042 NA
ENSG00000198517 MAFK 7975 -0.3776143 0.0012463155 0.8349042 NA
ENSG00000205763 RP9P NA 0.3568077 0.0012621256 0.8349042 NA
ENSG00000179163 FUCA1 2517 0.3866835 0.0012743644 0.8349042 NA
ENSG00000249715 FER1L5 90342 -0.9336025 0.0012876028 0.8349042 NA
$llBasalVs4
Low-responders Baseline vs 4 sets
ENSEMBL ID Gene name ENTREZID Log fold-change P-value Adjusted P-value labels
ENSG00000103507 BCKDK 10295 0.3577820 2.093069e-05 0.1491474 NA
ENSG00000143458 GABPB2 126626 -0.3943629 2.934603e-05 0.1491474 NA
ENSG00000105835 NAMPT 10135 0.4940640 3.452158e-05 0.1491474 NA
ENSG00000118298 CA14 23632 1.2067015 5.572125e-05 0.1491474 NA
ENSG00000170153 RNF150 57484 0.6365383 5.908350e-05 0.1491474 NA
ENSG00000177453 NIM1K 167359 1.6642045 6.037944e-05 0.1491474 NA
ENSG00000146278 PNRC1 10957 -0.4970777 6.282915e-05 0.1491474 NA
ENSG00000229931 ATXN1-AS1 NA -1.1191546 1.252427e-04 0.2432084 NA
ENSG00000067225 PKM 5315 0.4491092 1.317251e-04 0.2432084 NA
ENSG00000162433 AK4 205 0.9194336 1.875309e-04 0.2724934 NA
ENSG00000188825 LINC00910 100130581 -0.8145760 1.931663e-04 0.2724934 NA
ENSG00000165678 GHITM 27069 0.1945457 1.967817e-04 0.2724934 NA
ENSG00000165644 COMTD1 118881 0.4574684 2.207379e-04 0.2821539 NA
ENSG00000171747 LGALS4 3960 -0.7885530 3.095262e-04 0.3519180 NA
ENSG00000111669 TPI1 7167 0.3864333 3.517015e-04 0.3519180 NA
ENSG00000163933 RFT1 91869 -0.4356126 3.522775e-04 0.3519180 NA
ENSG00000067704 IARS2 55699 0.1919034 3.600292e-04 0.3519180 NA
ENSG00000141526 SLC16A3 9123 1.3022422 4.403330e-04 0.3635943 NA
ENSG00000144909 OSBPL11 114885 0.3925082 4.754764e-04 0.3635943 NA
ENSG00000117475 BLZF1 8548 -0.2669085 4.813656e-04 0.3635943 NA
$ll1Vs4
Low-responders 1 vs 4 sets
ENSEMBL ID Gene name ENTREZID Log fold-change P-value Adjusted P-value labels
ENSG00000101000 PROCR 10544 -1.0983144 4.715725e-05 0.7836120 NA
ENSG00000259868 AL163952.1 102723670 1.5549997 2.359479e-04 0.9996056 NA
ENSG00000198331 HYLS1 219844 0.6407774 8.839937e-04 0.9996056 NA
ENSG00000188825 LINC00910 100130581 -0.7195972 9.817186e-04 0.9996056 NA
ENSG00000239715 AC093627.2 105375115 1.4925838 1.052471e-03 0.9996056 NA
ENSG00000196388 INCA1 388324 0.3734889 1.181050e-03 0.9996056 NA
ENSG00000261934 PCDHGA9 56107 -1.0084618 1.227557e-03 0.9996056 NA
ENSG00000106804 C5 727 -1.3771154 1.255379e-03 0.9996056 NA
ENSG00000270673 YTHDF3-AS1 NA 1.1830918 1.274698e-03 0.9996056 NA
ENSG00000182575 NXPH3 11248 -0.6117099 1.441875e-03 0.9996056 NA
ENSG00000235257 ITGA9-AS1 101928153 -1.4820272 1.458171e-03 0.9996056 NA
ENSG00000134910 STT3A 3703 -0.2174099 1.515457e-03 0.9996056 NA
ENSG00000227467 LINC01537 101928555 -1.0544710 1.721499e-03 0.9996056 NA
ENSG00000090776 EFNB1 1947 -0.3662268 1.759498e-03 0.9996056 NA
ENSG00000129235 TXNDC17 84817 0.4317693 1.810457e-03 0.9996056 NA
ENSG00000275367 AC092111.1 NA 0.9985458 2.108988e-03 0.9996056 NA
ENSG00000113456 RAD1 5810 0.1814723 2.388151e-03 0.9996056 NA
ENSG00000115361 ACADL 33 0.4338935 2.418651e-03 0.9996056 NA
ENSG00000139908 TSSK4 283629 -0.6043970 2.426706e-03 0.9996056 NA
ENSG00000176390 CRLF3 51379 0.3147388 2.467640e-03 0.9996056 NA
$lhBasalVs1
Medium-responders Baseline vs 1 set
ENSEMBL ID Gene name ENTREZID Log fold-change P-value Adjusted P-value labels
ENSG00000145604 SKP2 6502 0.7167148 3.363424e-07 0.005589001 SKP2
ENSG00000049246 PER3 8863 -1.1928102 2.215261e-06 0.014317641 PER3
ENSG00000120733 KDM3B 51780 -0.3025514 2.584878e-06 0.014317641 KDM3B
ENSG00000117859 OSBPL9 114883 -0.3861375 3.867208e-06 0.014594777 OSBPL9
ENSG00000143889 HNRNPLL 92906 -0.3310508 4.391520e-06 0.014594777 HNRNPLL
ENSG00000171914 TLN2 83660 -0.5079292 5.994352e-06 0.016009503 TLN2
ENSG00000092096 SLC22A17 51310 0.8106111 6.785462e-06 0.016009503 SLC22A17
ENSG00000172932 ANKRD13D 338692 0.4946082 8.198588e-06 0.016009503 ANKRD13D
ENSG00000104660 LEPROTL1 23484 0.4629737 9.552271e-06 0.016009503 LEPROTL1
ENSG00000129521 EGLN3 112399 0.8049995 1.014878e-05 0.016009503 EGLN3
ENSG00000178409 BEND3 57673 0.7133636 1.059785e-05 0.016009503 BEND3
ENSG00000258461 AC012651.1 NA -0.4406476 1.250624e-05 0.017318013 AC012651.1
ENSG00000118689 FOXO3 2309 -0.8653564 1.508577e-05 0.017892765 FOXO3
ENSG00000121671 CRY2 1408 -0.3282430 1.589111e-05 0.017892765 CRY2
ENSG00000132326 PER2 8864 -1.2704906 1.615162e-05 0.017892765 PER2
ENSG00000071282 LMCD1 29995 -0.5818202 1.851404e-05 0.019227983 LMCD1
ENSG00000065970 FOXJ2 55810 -0.3363774 2.115112e-05 0.020674599 FOXJ2
ENSG00000124177 CHD6 84181 -0.4299687 2.362769e-05 0.021024104 CHD6
ENSG00000148339 SLC25A25 114789 -1.2075493 2.403911e-05 0.021024104 SLC25A25
ENSG00000173320 STOX2 56977 1.0603356 2.666871e-05 0.022157695 STOX2
$lhBasalVs4
Medium-responders Baseline vs 4 sets
ENSEMBL ID Gene name ENTREZID Log fold-change P-value Adjusted P-value labels
ENSG00000100603 SNW1 22938 -0.2616062 7.507284e-07 0.01187811 SNW1
ENSG00000196878 LAMB3 3914 1.4365939 1.882963e-06 0.01187811 LAMB3
ENSG00000121671 CRY2 1408 -0.3651942 2.144450e-06 0.01187811 CRY2
ENSG00000132326 PER2 8864 -1.4038004 3.303004e-06 0.01372150 PER2
ENSG00000215158 AC138409.2 NA -1.0104764 8.655706e-06 0.02496774 AC138409.2
ENSG00000134297 PLEKHA8P1 NA 0.9715068 9.015250e-06 0.02496774 PLEKHA8P1
ENSG00000175155 YPEL2 388403 0.5001611 1.068949e-05 0.02537533 YPEL2
ENSG00000020256 ZFP64 55734 0.4893317 1.394027e-05 0.02613278 ZFP64
ENSG00000118689 FOXO3 2309 -0.8666200 1.490782e-05 0.02613278 FOXO3
ENSG00000171621 SPSB1 80176 -1.0137618 1.572653e-05 0.02613278 SPSB1
ENSG00000049246 PER3 8863 -1.0567671 2.002753e-05 0.03025432 PER3
ENSG00000198218 QRICH1 54870 -0.1934987 2.434858e-05 0.03160221 QRICH1
ENSG00000204444 APOM 55937 -0.5630813 2.507725e-05 0.03160221 APOM
ENSG00000143393 PI4KB 5298 -0.2050077 2.746954e-05 0.03160221 PI4KB
ENSG00000164089 ETNPPL 64850 1.6694455 3.126509e-05 0.03160221 ETNPPL
ENSG00000145604 SKP2 6502 0.5745901 3.269894e-05 0.03160221 SKP2
ENSG00000128016 ZFP36 7538 -1.1468181 3.271954e-05 0.03160221 ZFP36
ENSG00000065970 FOXJ2 55810 -0.3273154 3.423240e-05 0.03160221 FOXJ2
ENSG00000171914 TLN2 83660 -0.4561842 3.892539e-05 0.03404333 TLN2
ENSG00000221817 PPP3CB-AS1 101929145 -0.7210390 4.307416e-05 0.03578817 PPP3CB-AS1
$lh1Vs4
Medium-responders 1 vs 4 sets
ENSEMBL ID Gene name ENTREZID Log fold-change P-value Adjusted P-value labels
ENSG00000218073 AL021407.3 NA 2.0350074 1.437948e-05 0.2296792 NA
ENSG00000049860 HEXB 3074 -0.3294538 2.764389e-05 0.2296792 NA
ENSG00000088543 C3orf18 51161 0.4040884 8.684265e-05 0.4252517 NA
ENSG00000068697 LAPTM4A 9741 -0.1840266 1.065912e-04 0.4252517 NA
ENSG00000260296 AC095057.3 NA 1.2960755 1.546834e-04 0.4252517 NA
ENSG00000166482 MFAP4 4239 -0.6521111 1.754570e-04 0.4252517 NA
ENSG00000248996 AC145098.1 NA 1.5102374 2.177369e-04 0.4252517 NA
ENSG00000269956 MKNK1-AS1 100507423 1.5191992 2.242062e-04 0.4252517 NA
ENSG00000166788 SAAL1 113174 0.4520025 2.797067e-04 0.4252517 NA
ENSG00000178175 ZNF366 167465 0.7107727 2.807466e-04 0.4252517 NA
ENSG00000121903 ZSCAN20 7579 0.7530707 3.070421e-04 0.4252517 NA
ENSG00000148344 PTGES 9536 -1.4466517 3.602023e-04 0.4252517 NA
ENSG00000087086 FTL 2512 -0.3405158 3.785628e-04 0.4252517 NA
ENSG00000272669 AL021707.6 NA 1.1976173 3.982583e-04 0.4252517 NA
ENSG00000226332 AL354836.1 NA -0.7992321 4.023136e-04 0.4252517 NA
ENSG00000112137 PHACTR1 221692 0.6768344 4.499209e-04 0.4252517 NA
ENSG00000184619 KRBA2 124751 1.0435474 4.538703e-04 0.4252517 NA
ENSG00000144120 TMEM177 80775 0.3530154 4.843389e-04 0.4252517 NA
ENSG00000260686 AC008669.1 NA 1.1465068 4.862359e-04 0.4252517 NA
ENSG00000253485 PCDHGA5 56110 -1.0960252 5.233906e-04 0.4348591 NA
$hhBasalVs1
High-responders Baseline vs 1 set
ENSEMBL ID Gene name ENTREZID Log fold-change P-value Adjusted P-value labels
ENSG00000163738 MTHFD2L 441024 -0.4346979 4.080866e-05 0.5071873 NA
ENSG00000198133 TMEM229B 161145 1.3639728 7.052866e-05 0.5071873 NA
ENSG00000107485 GATA3 2625 1.0671820 1.080818e-04 0.5071873 NA
ENSG00000203808 BVES-AS1 154442 -0.3257393 1.220888e-04 0.5071873 NA
ENSG00000124574 ABCC10 89845 0.2927855 2.768761e-04 0.5652657 NA
ENSG00000170142 UBE2E1 7324 -0.1812900 2.999777e-04 0.5652657 NA
ENSG00000269973 AC010969.2 NA 1.2074584 3.325905e-04 0.5652657 NA
ENSG00000198382 UVRAG 7405 -0.2498498 3.774890e-04 0.5652657 NA
ENSG00000138279 ANXA7 310 -0.1704141 3.975185e-04 0.5652657 NA
ENSG00000237298 TTN-AS1 100506866 -0.2890136 4.130979e-04 0.5652657 NA
ENSG00000241506 PSMC1P1 151645 -0.3507469 4.420839e-04 0.5652657 NA
ENSG00000102978 POLR2C 5432 -0.2177528 4.708028e-04 0.5652657 NA
ENSG00000012171 SEMA3B 7869 0.4571562 4.904573e-04 0.5652657 NA
ENSG00000196123 KIAA0895L 653319 0.4155749 5.096112e-04 0.5652657 NA
ENSG00000159640 ACE 1636 0.5764126 5.500382e-04 0.5652657 NA
ENSG00000105137 SYDE1 85360 0.4803123 5.937038e-04 0.5652657 NA
ENSG00000166341 DCHS1 8642 0.4790866 6.354315e-04 0.5652657 NA
ENSG00000118298 CA14 23632 1.0277066 7.389403e-04 0.5652657 NA
ENSG00000260852 FBXL19-AS1 283932 0.7534122 7.660094e-04 0.5652657 NA
ENSG00000155760 FZD7 8324 -0.4935716 8.075141e-04 0.5652657 NA
$hhBasalVs4
High-responders Baseline vs 4 sets
ENSEMBL ID Gene name ENTREZID Log fold-change P-value Adjusted P-value labels
ENSG00000135924 DNAJB2 3300 -0.3225192 4.574562e-07 0.006379536 DNAJB2
ENSG00000077009 NMRK2 27231 1.2852016 7.678324e-07 0.006379536 NMRK2
ENSG00000041880 PARP3 10039 0.4069532 2.286676e-06 0.012665899 PARP3
ENSG00000118298 CA14 23632 1.4213099 4.599813e-06 0.017499001 CA14
ENSG00000187699 C2orf88 84281 0.5899930 5.363784e-06 0.017499001 C2orf88
ENSG00000213199 ASIC3 9311 0.5693084 6.318469e-06 0.017499001 ASIC3
ENSG00000156265 MAP3K7CL 56911 -0.8894052 1.102713e-05 0.023126324 MAP3K7CL
ENSG00000117450 PRDX1 5052 -0.3274787 1.113381e-05 0.023126324 PRDX1
ENSG00000135537 AFG1L 246269 0.6685325 1.605258e-05 0.028161215 AFG1L
ENSG00000159352 PSMD4 5710 -0.2338607 1.696446e-05 0.028161215 PSMD4
ENSG00000138279 ANXA7 310 -0.2057951 2.624342e-05 0.028161215 ANXA7
ENSG00000104852 SNRNP70 6625 -0.3512945 2.854709e-05 0.028161215 SNRNP70
ENSG00000173113 TRMT112 51504 -0.3998957 2.996597e-05 0.028161215 TRMT112
ENSG00000180867 PDIA3P1 NA 0.7975326 3.257907e-05 0.028161215 PDIA3P1
ENSG00000136826 KLF4 9314 -0.8579446 3.424681e-05 0.028161215 KLF4
ENSG00000234281 LANCL1-AS1 NA -0.7800724 3.541456e-05 0.028161215 LANCL1-AS1
ENSG00000187498 COL4A1 1282 0.8247442 3.573222e-05 0.028161215 COL4A1
ENSG00000251602 AL928654.1 100507437 0.7167614 3.695756e-05 0.028161215 AL928654.1
ENSG00000067225 PKM 5315 0.4886361 3.715772e-05 0.028161215 PKM
ENSG00000154553 PDLIM3 27295 -0.5539931 3.880111e-05 0.028161215 PDLIM3
$hh1Vs4
High-responders 1 vs 4 sets
ENSEMBL ID Gene name ENTREZID Log fold-change P-value Adjusted P-value labels
ENSG00000100234 TIMP3 7078 -0.6618710 4.455894e-05 0.6248607 NA
ENSG00000027644 INSRR 3645 -1.7977718 1.665942e-04 0.6248607 NA
ENSG00000154134 ROBO3 64221 -0.7603679 1.678863e-04 0.6248607 NA
ENSG00000122873 CISD1 55847 0.2779992 2.406105e-04 0.6248607 NA
ENSG00000271855 AC073195.1 NA -1.6429918 2.704983e-04 0.6248607 NA
ENSG00000222112 RN7SKP16 106480837 1.2782063 3.007082e-04 0.6248607 NA
ENSG00000152475 ZNF837 116412 -1.4656571 3.131693e-04 0.6248607 NA
ENSG00000110934 BIN2 51411 1.5348528 3.135429e-04 0.6248607 NA
ENSG00000166482 MFAP4 4239 -0.5987242 4.263965e-04 0.6248607 NA
ENSG00000139874 SSTR1 6751 -1.2661438 4.623963e-04 0.6248607 NA
ENSG00000213463 SYNJ2BP 55333 0.2251943 4.970576e-04 0.6248607 NA
ENSG00000183401 CCDC159 126075 -0.4916470 4.988519e-04 0.6248607 NA
ENSG00000196378 ZNF34 80778 -0.4585302 5.267359e-04 0.6248607 NA
ENSG00000041880 PARP3 10039 0.2852183 5.311601e-04 0.6248607 NA
ENSG00000186310 NAP1L3 4675 -1.2688012 5.897903e-04 0.6248607 NA
ENSG00000154864 PIEZO2 63895 -1.9417950 6.113453e-04 0.6248607 NA
ENSG00000104738 MCM4 4173 0.2720460 7.214167e-04 0.6248607 NA
ENSG00000163738 MTHFD2L 441024 0.3486109 8.395199e-04 0.6248607 NA
ENSG00000151687 ANKAR 150709 0.4998526 8.673863e-04 0.6248607 NA
ENSG00000135537 AFG1L 246269 0.4875207 9.163581e-04 0.6248607 NA

Graphics RNAseq

PCA plot

# Prepare x_pca data.frame
x_pca <- x$counts |> 
  as.data.frame() |> 
  rownames_to_column("Genes") |> 
  pivot_longer(cols = -Genes, names_to = "ID", values_to = "values") |> 
  pivot_wider(id_cols = ID, names_from = Genes, values_from = values) |> 
  column_to_rownames("ID")

# Run PCA
pca <- prcomp(x_pca, scale. = T, center = T)

# Raw PCA from factoextra package
fviz_pca_ind(pca)

# Create PCA data.frame
pca_df <- as.data.frame(pca$x) |> 
  rownames_to_column("ID") |>
  left_join(as.data.frame(x$samples) |> rownames_to_column("ID"),, join_by("ID")) |> 
  mutate(Groups = case_when(group == "L1_L4" ~ "Low-responders",
                            group == "L1_H4" ~ "Medium-responders",
                            group == "H1_H4" ~ "High-responders")) |> 
  mutate(Groups = factor(Groups, levels = c("Low-responders", "Medium-responders", "High-responders")))

# Text configuration for graph
config_text <- element_text(size = 14, colour = "black", family = "serif")

pca_df |>
  ggplot(aes(x = PC1, y = PC2, fill = Groups)) +
  geom_point(shape = 21, size = 3, show.legend = T, alpha = 0.8) +
  theme_minimal() +
  labs(x = "PC1 (28.1%)", y = "PC2 (6.7%)", fill = "") + 
  geom_hline(yintercept = 0, linetype = 2) +
  geom_vline(xintercept = 0, linetype = 2) +
  theme(text = config_text,
        axis.text = config_text,
        legend.text = config_text)

Volcano function

fun_volcano <- function(data, title){
  
  volcano = EnhancedVolcano(data, 
                            lab = data$external_gene_name,
                            x = "logFC",
                            y = "adj.P.Val",
                            pCutoff = 0.05,
                            FCcutoff = 1,
                            # pCutoffCol = "adj.P.Val",
                            ylab = bquote(~-Log[10] ~ Adjusted~italic(P)-Values),
                            title = title,
                            subtitle = NULL,
                            axisLabSize = 12,
                            titleLabSize = 12,
                            labSize = 3,
                            caption = NULL,
                            legendPosition = "NA",
                            max.overlaps = Inf,
                            drawConnectors = T,
                            lengthConnectors = unit(0.05, "npc"),
                            arrowheads = F,
                            ylim = c(0, 5))
  
  return(volcano)
}

Volcano plots

# Create names
groups_names <- paste(rep(c("Low-responders", "Medium-responders", "High-responders"), each = 3), c("Baseline vs 1 set", "Baseline vs 4 sets", "1 vs 4 sets")) 
# map volcano function
volcano_plots <- map2(list(llBasalVs1, llBasalVs4, ll1Vs4, lhBasalVs1, lhBasalVs4, lh1Vs4, hhBasalVs1, hhBasalVs4, hh1Vs4), groups_names, \(x, y) fun_volcano(x, y))

cowplot:::plot_grid(plotlist = volcano_plots,
                    labels = LETTERS[1:9],
                    nrow = 3, 
                    ncol = 3)

GO annotation from clusterprofiler

Function GO over-representation annotation

fun_go <- function(data, ...){
  
  df <- DE_list_filtered[[data]] |> 
    filter(!is.na(ENTREZID)) |> 
    filter(`Adjusted P-value` <= 0.05)
  
  gene <- df$ENTREZID
  
  if(length(gene) >= 1){
  
  go <- enrichGO(gene = gene,
                 universe = DE_list_filtered[[data]][["ENTREZID"]],
                 OrgDb = org.Hs.eg.db,
                 ont = "BP",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.05,
                 qvalueCutoff = 0.2,
                 readable = T)
  
  go1 <- go |> 
    as_tibble() |> 
    rename("Gene ratio" = GeneRatio, "BG ratio" = BgRatio, "P-value" = pvalue, "Adjusted P-value" = p.adjust, "Q-value" = qvalue, "Gene ID" = geneID)
  
  return(go1)
  
  } else {
    
    return("No genes found")
    
  }
  
}

Go-over represetation annotation

# Create names
go_names <- list("llBasalVs1", "llBasalVs4", "ll1Vs4", "lhBasalVs1", "lhBasalVs4", "lh1Vs4", "hhBasalVs1", "hhBasalVs4", "hh1Vs4")

# map GO function
go_enrichment_list <- map(go_names, \(x) fun_go(x)) |> 
  set_names(go_names)

go_enrichment_list |>
  map(\(x) as.data.frame(x)) |> 
  map(\(x) head(x, 20)) |> 
  map2(groups_names, \(x, header) gt(x) |> tab_header(title = header))
$llBasalVs1
Low-responders Baseline vs 1 set
x
No genes found
$llBasalVs4
Low-responders Baseline vs 4 sets
x
No genes found
$ll1Vs4
Low-responders 1 vs 4 sets
x
No genes found
$lhBasalVs1
Medium-responders Baseline vs 1 set
ID Description Gene ratio BG ratio RichFactor FoldEnrichment zScore P-value Adjusted P-value Q-value Gene ID Count
GO:0048511 rhythmic process 11/76 233/12312 0.04721030 7.648069 8.074055 1.709722e-07 0.0002542357 0.0002426005 PER3/FOXO3/CRY2/PER2/PER1/NR1D2/NFIL3/CIART/BHLHE41/TEF/SGPL1 11
GO:0032922 circadian regulation of gene expression 6/76 61/12312 0.09836066 15.934426 9.215129 1.920561e-06 0.0014279372 0.0013625876 PER3/CRY2/PER2/PER1/CIART/BHLHE41 6
GO:0007623 circadian rhythm 8/76 160/12312 0.05000000 8.100000 7.124070 6.084936e-06 0.0026807268 0.0025580432 PER3/CRY2/PER2/PER1/NR1D2/NFIL3/CIART/BHLHE41 8
GO:0043153 entrainment of circadian clock by photoperiod 4/76 21/12312 0.19047619 30.857143 10.791902 7.409589e-06 0.0026807268 0.0025580432 PER3/CRY2/PER2/PER1 4
GO:0009648 photoperiodism 4/76 22/12312 0.18181818 29.454545 10.527392 9.013876e-06 0.0026807268 0.0025580432 PER3/CRY2/PER2/PER1 4
GO:0009649 entrainment of circadian clock 4/76 24/12312 0.16666667 27.000000 10.047826 1.297187e-05 0.0032148611 0.0030677328 PER3/CRY2/PER2/PER1 4
GO:0042921 nuclear receptor-mediated glucocorticoid signaling pathway 3/76 11/12312 0.27272727 44.181818 11.291738 3.599440e-05 0.0076462390 0.0072963085 CRY2/PER1/NEDD4 3
GO:0031958 nuclear receptor-mediated corticosteroid signaling pathway 3/76 12/12312 0.25000000 40.500000 10.788694 4.777946e-05 0.0088810080 0.0084745682 CRY2/PER1/NEDD4 3
GO:0071383 cellular response to steroid hormone stimulus 7/76 170/12312 0.04117647 6.670588 5.867339 8.299393e-05 0.0137124422 0.0130848916 SKP2/FOXO3/CRY2/PER1/NEDD4/NR3C2/ZFP36 7
GO:0042752 regulation of circadian rhythm 5/76 89/12312 0.05617978 9.101124 6.044844 2.174350e-04 0.0323325905 0.0308528879 PER3/CRY2/PER2/PER1/NR1D2 5
$lhBasalVs4
Medium-responders Baseline vs 4 sets
ID Description Gene ratio BG ratio RichFactor FoldEnrichment zScore P-value Adjusted P-value Q-value Gene ID Count
GO:0042754 negative regulation of circadian rhythm 3/25 10/12312 0.30000000 147.74400 20.939496 8.792305e-07 0.0007227275 0.0006062063 CRY2/PER2/SFPQ 3
GO:0043153 entrainment of circadian clock by photoperiod 3/25 21/12312 0.14285714 70.35429 14.347726 9.602116e-06 0.0029908237 0.0025086305 CRY2/PER2/PER3 3
GO:0009648 photoperiodism 3/25 22/12312 0.13636364 67.15636 14.008795 1.110334e-05 0.0029908237 0.0025086305 CRY2/PER2/PER3 3
GO:0009649 entrainment of circadian clock 3/25 24/12312 0.12500000 61.56000 13.395059 1.455389e-05 0.0029908237 0.0025086305 CRY2/PER2/PER3 3
GO:0042752 regulation of circadian rhythm 4/25 89/12312 0.04494382 22.13393 9.025685 2.873656e-05 0.0047242909 0.0039626208 CRY2/PER2/PER3/SFPQ 4
GO:0045646 regulation of erythrocyte differentiation 3/25 42/12312 0.07142857 35.17714 10.007644 8.058091e-05 0.0106463988 0.0089299413 FOXO3/ZFP36/LDB1 3
GO:0048511 rhythmic process 5/25 233/12312 0.02145923 10.56824 6.651021 9.066276e-05 0.0106463988 0.0089299413 CRY2/PER2/FOXO3/PER3/SFPQ 5
GO:0032922 circadian regulation of gene expression 3/25 61/12312 0.04918033 24.22033 8.200517 2.462733e-04 0.0253045765 0.0212248657 CRY2/PER2/PER3 3
GO:0007623 circadian rhythm 4/25 160/12312 0.02500000 12.31200 6.496348 2.808376e-04 0.0256498344 0.0215144597 CRY2/PER2/PER3/SFPQ 4
GO:0071383 cellular response to steroid hormone stimulus 4/25 170/12312 0.02352941 11.58776 6.270142 3.538365e-04 0.0290853631 0.0243960978 CRY2/FOXO3/SKP2/ZFP36 4
$lh1Vs4
Medium-responders 1 vs 4 sets
x
No genes found
$hhBasalVs1
High-responders Baseline vs 1 set
x
No genes found
$hhBasalVs4
High-responders Baseline vs 4 sets
ID Description Gene ratio BG ratio RichFactor FoldEnrichment zScore P-value Adjusted P-value Q-value Gene ID Count
GO:0042026 protein refolding 3/32 23/12312 0.1304348 50.18478 12.05199 2.726577e-05 0.02344856 0.02143951 DNAJB2/SNRNP70/HSPB1 3
$hh1Vs4
High-responders 1 vs 4 sets
x
No genes found

Plier analysis

Data preparation

z <- cpm(x, log=TRUE)

# Calculate the means of each row in the dataset
z.means <- rowMeans(z)

# Select only genes that have a correspondence in the z dataset
gene <- gene.list[rownames(z), ]

# Applies a function to each item in 1:nrow(z) to find the location of the maximum value.
ensembl_to_gene <- tapply(1:nrow(z), as.character(gene$external_gene_name), function(i){
    i[which.max(z.means[i])]})

# Dataset that I'll be working with.
y <- z[ensembl_to_gene, ]

rownames(y) <- names(ensembl_to_gene)

Load Plier pathways

# Load dataset
load("../0_data/MetaMEx_matrix.Rds")
data(bloodCellMarkersIRISDMAP)
data(canonicalPathways)

allPaths <- combinePaths(bloodCellMarkersIRISDMAP, canonicalPathways, MetaMEx_matrix)

Prepare gene and pathway sets

# Find genes in the pathways that are in our filtered data set
cm.genes <- commonRows(allPaths, y)

length(cm.genes)
[1] 7221
# Normalize expression using z-score
yN <- PLIER::rowNorm(y)
yN <- data.frame(yN)

# Compute principal components
PLIER::num.pc(yN) 
[1] 39

Load Plier results

  • By loading this, you don’t need to re-run the code below.
load("./0_data/Plier_results_object_trimmed_final.Rdata")

Run PLIER

# Singular value decomposition of a matrix (SVD) decomposition on the same set of genes.
# Additional information about SVD can be found here: https://web.mit.edu/be.400/www/SVD/Singular_Value_Decomposition.htm
yN.svd <- svd(yN[cm.genes,])

# Pathway matrix pseudo-inverse . This is used to pre-select pathways for optimization.Chat 

Chat = computeChat(allPaths[cm.genes,])

plierResults <- PLIER(as.matrix(yN[cm.genes, ]), allPaths[cm.genes, ], k=39, svdres=yN.svd, Chat = Chat)

Plot matrix

dim(yN)
[1] 16406    81
plotU(plierResults, auc.cutoff = 0.70, fdr.cutoff = 0.1, top = 3)

Interpret PLIER results

#create a data.frame with the Plier results
scores <- as.data.frame(t(plierResults$B))

for (i in 1:length(colnames(scores))){
  colnames(scores)[i] <- paste(c("LV",i), collapse = "")
}

my.design <- design.data 

scores$Group <- factor(my.design$Group[match(rownames(scores), rownames(my.design))])
scores$TimePoint <- factor(my.design$TimePoint[match(rownames(scores), rownames(my.design))])
scores$Set <- factor(my.design$Set[match(rownames(scores), rownames(my.design))])
scores$Sex <- factor(my.design$Sex[match(rownames(scores), rownames(my.design))])
scores$subjectID <- rownames(scores)

# Data wrangling for the score dataset
scores_final <- scores |> 
  mutate(Group = case_when(Group == "L1_L4" ~ "LL",
                           Group == "L1_H4" ~ "LH",
                           Group == "H1_H4" ~ "HH")) |> 
  mutate(Combined = paste(TimePoint, Set, sep = "_"),
         Combined = case_when(Combined == "pre_1" ~ "Baseline",
                              Combined == "post_1" ~ "Post_1set",
                              Combined == "post_4" ~ "Post_4set")) |> 
  mutate(ID = as.numeric(map(str_extract_all(subjectID, "[:number:]"), \(x) str_c(x, collapse = ""))))

# Run mixed model analyses for all LVs (n = 39)
model_RNA <- names(scores_final)[1:39] |> 
  map(\(x) lmer(formula = as.formula(paste0(x, " ~ Group * Combined + (1 | ID)")), data = scores_final)) |> map(\(x) car::Anova(x, test.statistic = "F")) |> 
  map_df(\(x) x) |> 
  mutate(LV = rep(seq(1, 39), each = 3), Effect = rep(c("Group", "Set", "Group:Set"), 39))

# Filter LVs with group:set interaction
model_rna_df <- model_RNA |> 
  filter(Effect == "Group:Set" & `Pr(>F)` <= 0.05)

# Select LVs of interest
LVs <- paste0("LV", model_rna_df[["LV"]])

# Final dataset
RNAseq_df <- scores_final |> 
  dplyr::select(ID, Group, Combined, all_of(LVs)) |> 
  dplyr::rename("Set" = Combined) |> 
  mutate(Set = case_when(Set == "Baseline" ~ "Baseline",
                         Set == "Post_1set" ~ "1 Set",
                         Set == "Post_4set" ~ "4 Sets")) |> 
  mutate(Set = factor(Set, levels = c("Baseline", "1 Set", "4 Sets"))) |> 
  mutate(Group = factor(Group, levels = c("LL", "LH", "HH"), labels = c("Low-responders", "Medium-responders", "High-responders"))) |> 
  arrange(Group, Set)

gt(RNAseq_df)
ID Group Set LV3 LV8 LV11 LV13 LV27 LV32 LV38
102 Low-responders Baseline -0.461586424 0.8894738481 -0.12340442 0.508480173 0.485320048 -1.803526449 1.129070872
12 Low-responders Baseline 0.359095588 0.8882311941 0.33872514 -0.311425929 0.140069166 -0.597810030 0.577088450
30 Low-responders Baseline 0.165720424 -0.1516787365 -0.64836282 -0.451114127 -0.097700458 0.698137599 -0.624720500
35 Low-responders Baseline 0.091675528 -0.1812164228 -0.07827111 0.676322276 -0.485598779 -0.881092536 -0.110975882
59 Low-responders Baseline 0.059861689 0.5157809574 0.34634407 -0.715988679 -0.821810849 0.085191649 -0.556102877
5 Low-responders Baseline -0.084349255 1.4365378719 -0.06770246 0.472668845 -0.263039595 -0.306187762 0.052146211
77 Low-responders Baseline 0.341389383 -0.3460115321 0.26911637 -0.658588247 0.436643548 0.471478244 0.542628876
90 Low-responders Baseline -1.839555113 0.0861968728 -0.05064288 0.800528890 0.489674220 -0.978101820 0.399021254
94 Low-responders Baseline 0.154525953 -0.9813190551 -0.14411191 -0.468985595 -0.342966561 1.637656668 -0.789799833
102 Low-responders 1 Set -0.355749956 0.3783556484 -0.05277054 -0.670276750 1.307757652 0.540977271 1.481922672
30 Low-responders 1 Set 0.069628545 1.3037979779 -0.35211836 -1.056288033 0.243141867 0.498964363 -0.781418998
59 Low-responders 1 Set 0.623139456 0.5078404653 0.08489580 -0.232708278 -1.063181253 0.109547849 -0.525626919
94 Low-responders 1 Set -0.116009196 -0.8716610766 0.33510814 -0.936870912 0.017316455 1.075728030 -0.706625719
12 Low-responders 1 Set -0.284969166 0.9983143568 0.39339712 0.028799849 0.925363732 0.589205117 0.449727953
35 Low-responders 1 Set 0.177057554 0.3156103545 0.50471377 -0.909538809 -0.424256360 0.203638486 -0.822785609
5 Low-responders 1 Set -0.249147578 1.1972058760 -0.01506670 0.882872605 -0.045161807 0.809628404 0.140832563
77 Low-responders 1 Set 0.222860502 0.2807024635 0.47538791 0.130414440 -0.228579025 -1.539009595 0.364324347
90 Low-responders 1 Set -2.032719876 -0.0279682824 0.48612397 0.200966255 0.756753051 0.657799412 0.213932120
12 Low-responders 4 Sets -0.113579706 0.1986182506 0.47513523 0.180847534 0.669104311 0.461472055 0.381277373
35 Low-responders 4 Sets -0.002175232 1.0184304475 0.60066684 -0.847959717 0.175200025 0.124369513 -0.583682035
5 Low-responders 4 Sets -0.209507283 0.8358493909 -0.01249616 0.634341946 -0.007743541 0.736212018 0.266639966
77 Low-responders 4 Sets 0.164739566 -0.5821048975 0.39068933 0.504944606 0.275117843 -0.984675191 0.369947375
90 Low-responders 4 Sets -1.798410293 -0.4889108235 0.39252995 0.221988484 0.668234404 0.449405286 0.236223800
102 Low-responders 4 Sets -0.426532418 -0.0738016507 0.04444962 -0.699984626 1.020427413 0.382373580 1.556101955
30 Low-responders 4 Sets 0.063976860 0.8318398043 0.05490554 -1.081060079 -0.188920388 0.312272715 -0.847864458
59 Low-responders 4 Sets 0.314686135 0.1927773758 -0.17430274 -0.847636468 -0.559355329 0.174697983 -0.486411250
94 Low-responders 4 Sets 0.040510310 -1.4347184959 0.17813722 -0.633217982 -0.098115543 1.252927015 -0.840475081
105 Medium-responders Baseline 0.646056592 0.4475398493 -0.18797995 -0.636748521 -0.926053480 0.804247506 -0.006388977
106 Medium-responders Baseline -0.023722570 -0.0026040539 -0.29887412 -0.463387466 -0.281419557 0.776468996 -1.302322094
107 Medium-responders Baseline 0.012558347 0.8982902322 0.27593768 0.358132204 0.847796296 0.703719858 -0.138272806
110 Medium-responders Baseline -0.339829935 0.3831177235 1.06709574 -0.130689319 -0.417460635 0.471470571 -1.427321886
114 Medium-responders Baseline -0.089359403 -1.3441590778 0.69927292 -0.394676511 -0.209037700 -0.401582671 0.108422552
17 Medium-responders Baseline 0.149361864 -0.0572958781 0.04793061 -0.343134578 0.222260112 1.701452824 1.122141431
50 Medium-responders Baseline 0.200611740 1.0569184873 0.79007972 0.033106948 -0.451793832 0.003670448 -1.743938767
76 Medium-responders Baseline -0.085688691 -0.2758317279 -0.17164339 1.259564203 1.260831857 0.360165901 -0.071720842
98 Medium-responders Baseline -0.059761290 0.2071664873 1.36427556 -0.574674140 0.520655394 -0.517820154 0.261049209
110 Medium-responders 1 Set 0.496517896 1.2571270832 0.18095727 1.096014967 -1.074713461 -1.414088799 -0.512240950
17 Medium-responders 1 Set 0.267160600 0.3936156168 -0.06171895 0.142987533 0.232885980 -0.199206779 1.282730110
50 Medium-responders 1 Set 0.220017462 1.9591383736 0.13996344 -0.654385744 -0.617542453 -0.769095132 -1.495873877
105 Medium-responders 1 Set 0.499938427 0.6105145273 0.35149379 1.009082957 -1.332301922 -0.698669700 0.366804448
106 Medium-responders 1 Set 0.095879918 -0.0439601883 -0.74368377 0.649290638 -0.433553380 -0.717400349 -0.659109682
107 Medium-responders 1 Set 0.347908774 0.7235800489 -0.18250858 1.169369917 -0.249720512 -0.484312733 -0.113368636
114 Medium-responders 1 Set 0.333892614 -0.9356047044 0.83521411 -0.649124481 -1.251827737 -0.472585934 -0.265572767
76 Medium-responders 1 Set 0.050663627 -0.6827226613 0.11509006 1.557120518 0.805346107 -0.087900872 0.132383847
98 Medium-responders 1 Set 0.051444174 -0.5441520920 1.14647488 0.073669135 0.341739456 -1.120811361 0.759168448
105 Medium-responders 4 Sets 0.250070058 0.1832833847 0.35325979 1.096272997 -1.115740894 -0.821048026 0.253627985
106 Medium-responders 4 Sets -0.125421553 -1.5200570075 -0.50352115 1.144852984 -0.032347013 -0.516016921 -0.985606128
107 Medium-responders 4 Sets 0.244461357 -0.3943102615 -0.26619416 1.124123925 0.542397211 -0.383913525 0.244143465
114 Medium-responders 4 Sets 0.067168490 -1.4248300042 0.40004262 -0.255231232 -0.094924126 0.461040908 0.139896413
76 Medium-responders 4 Sets 0.099846224 -0.2324002429 0.12206759 1.598871394 0.734163412 -0.083385548 -0.099230016
98 Medium-responders 4 Sets -0.065193191 -0.5033268292 1.06318318 -0.187297333 0.352481361 -1.132673261 0.868002960
110 Medium-responders 4 Sets -0.037268857 -0.8221156564 0.26071988 1.343272740 -0.148379621 -0.753778061 -0.785448911
17 Medium-responders 4 Sets 0.228187387 -0.0078986586 0.05460092 0.138410996 -0.640094863 -0.722958634 0.850363544
50 Medium-responders 4 Sets 0.184186006 -0.0880514028 0.58749525 0.019885285 -0.747171423 -0.529335503 -1.811534090
104 High-responders Baseline -0.074165934 -0.1254202033 -0.16363153 -0.127286692 1.316957581 0.489128195 -0.210328824
27 High-responders Baseline -0.363829839 -0.1641205852 -0.03629102 -0.648760702 0.525216786 -0.105801086 0.967516667
32 High-responders Baseline -0.124125234 0.9291745792 -0.75312646 -0.006708657 0.062680676 0.457340692 0.678506269
36 High-responders Baseline 0.226418191 0.6595776235 -0.46104679 0.408279762 -0.639942721 -0.899241827 -0.405129107
38 High-responders Baseline -0.448632505 -0.1646702369 -0.76073344 0.088069861 -1.082709903 -0.107321945 0.344415507
6 High-responders Baseline 0.145756707 0.2966271813 0.27094260 -0.782164839 -0.514831499 0.429333332 -0.321768039
80 High-responders Baseline -0.013090050 -0.1733896194 -0.05766169 -0.359529010 0.267436417 1.030815680 0.735862394
86 High-responders Baseline 0.092836295 0.5394826390 -0.32826650 -0.251198200 -0.110230287 1.243968467 -0.802577100
97 High-responders Baseline 0.538528016 -0.8746712258 -0.54374769 -0.492779555 0.292325462 0.541919812 -1.269619286
36 High-responders 1 Set 0.561371025 -0.1777275554 -0.77945166 0.707793935 -0.273840138 -0.621511086 0.297911926
80 High-responders 1 Set -0.060459524 -0.5134181734 0.21799476 0.150500199 0.460712413 0.163848067 0.860058832
97 High-responders 1 Set 0.542356175 -1.0319792285 -0.70085250 -0.905983581 0.014730829 0.178746779 -0.742754729
104 High-responders 1 Set -0.150572016 0.1665457035 -0.68056037 -0.380748468 1.287643352 0.006281901 -0.018468542
27 High-responders 1 Set -0.196646689 0.0009305934 -0.13651850 -1.425118993 0.458200557 0.202435879 1.226680894
32 High-responders 1 Set 0.041840697 0.3786504403 -0.30550408 -0.605125949 -0.110350508 -0.214844523 0.448833690
38 High-responders 1 Set -0.184697794 -0.9991895548 -1.08097154 0.753680499 -1.318629309 -0.318455899 0.233795107
6 High-responders 1 Set 0.023010449 -0.7623429530 0.24694233 -0.112481184 -0.348274160 0.656560456 -0.519996726
86 High-responders 1 Set 0.143941058 -0.5643131754 -0.42735375 1.045434923 -0.375606347 -0.125786452 -0.303809515
104 High-responders 4 Sets 0.204660828 0.0829011188 -0.74909413 -0.402323980 1.729644632 0.221037895 0.079096658
27 High-responders 4 Sets -0.206984730 0.0041405529 -0.08913098 -1.101114320 0.624009577 0.348318388 1.117737814
32 High-responders 4 Sets -0.178016887 -0.4945517437 -0.56675078 -0.004745324 0.361871855 0.167464188 1.014145692
38 High-responders 4 Sets -0.113417798 0.2275778892 -1.33067673 0.522982260 -0.656033969 -0.413085408 0.688668927
6 High-responders 4 Sets 0.097259184 -0.5526066833 0.19402050 -0.193303928 -0.311025456 0.511290233 -0.279623543
86 High-responders 4 Sets 0.156976279 -1.0448961947 -0.63620616 0.981821866 -0.271899821 -0.280459855 0.032577616
36 High-responders 4 Sets 0.353890325 -0.8043527611 -0.76778329 0.420968309 -0.354771440 -0.824752098 0.360040977
80 High-responders 4 Sets -0.047042557 0.5146949083 0.04669738 -0.606503976 0.263827249 -0.064540512 1.101047135
97 High-responders 4 Sets 0.538574266 -1.3297968853 -0.71134483 -0.919865945 -0.117280649 -0.309622198 -0.768001303

Posthoc for significant LVs

# Run mixed model
model_rna_posthoc <- names(RNAseq_df)[4:10] |>
  map(\(x) lmer(formula = as.formula(paste(x, "~ Group * Set + (1 | ID)")), data = RNAseq_df)) |> 
  # map(\(x) car::Anova(x, test.statistic = "F")) |> 
  set_names(names(RNAseq_df)[4:10]) 

# Run posthoc mixed model
posthoc_rna <- map(1:7, \(x) emmeans(model_rna_posthoc[[x]], ~ Group * Set, adjust = "tukey")) |>
  (\(data) map(1:7, \(x) pairs(data[[x]], reverse = T) ))() |> 
  set_names(names(RNAseq_df)[4:10]) 

# Set confidence interval for the pairwise comparison
posthoc_rna_ci <- map(1:7, \(x) confint(posthoc_rna[[x]])) |> 
  set_names(names(RNAseq_df)[4:10])

# Create data frame with the results
posthoc_final <- posthoc_rna |> 
  map(\(x) as_tibble(x)) |> 
  (\(data) map(1:7, \(x) left_join(data[[x]], posthoc_rna_ci[[x]])))() |> 
  map(\(x) x |> mutate(significance = case_when(p.value <= 0.1 & p.value > 0.05 ~ "NS",
                                                p.value <= 0.05 ~ "**",
                                                p.value > 0.05 ~ "NS")) |> 
      mutate_at(.vars = c(2:8), \(y) round(y, 3)) |> 
      filter(p.value <= 0.05)) |>  
  set_names(names(RNAseq_df)[4:10])

# Filtered by p.values <= 0.05
posthoc_final |> 
  map(\(x) head(x, 20)) |> 
  imap(\(x, header) gt(x) |> tab_header(title = header))
$LV3
LV3
contrast estimate SE df t.ratio p.value lower.CL upper.CL significance
$LV8
LV8
contrast estimate SE df t.ratio p.value lower.CL upper.CL significance
(Medium-responders 4 Sets) - (Medium-responders 1 Set) -0.839 0.227 48 -3.693 0.015 -1.576 -0.101 **
$LV11
LV11
contrast estimate SE df t.ratio p.value lower.CL upper.CL significance
(High-responders Baseline) - (Medium-responders Baseline) -0.713 0.205 34.155 -3.486 0.033 -1.390 -0.036 **
(High-responders 1 Set) - (Medium-responders Baseline) -0.804 0.205 34.155 -3.927 0.010 -1.480 -0.127 **
(High-responders 4 Sets) - (Medium-responders Baseline) -0.911 0.205 34.155 -4.451 0.002 -1.588 -0.234 **
(High-responders 4 Sets) - (Low-responders 1 Set) -0.719 0.205 34.155 -3.513 0.030 -1.396 -0.042 **
(High-responders 4 Sets) - (Medium-responders 1 Set) -0.710 0.205 34.155 -3.471 0.034 -1.387 -0.033 **
(High-responders 4 Sets) - (Low-responders 4 Sets) -0.729 0.205 34.155 -3.562 0.027 -1.406 -0.052 **
(High-responders 4 Sets) - (Medium-responders 4 Sets) -0.742 0.205 34.155 -3.628 0.023 -1.419 -0.066 **
$LV13
LV13
contrast estimate SE df t.ratio p.value lower.CL upper.CL significance
(Medium-responders 4 Sets) - (Medium-responders Baseline) 0.768 0.207 48 3.719 0.014 0.097 1.439 **
$LV27
LV27
contrast estimate SE df t.ratio p.value lower.CL upper.CL significance
(Medium-responders 1 Set) - (Medium-responders Baseline) -0.461 0.129 48 -3.562 0.022 -0.881 -0.041 **
$LV32
LV32
contrast estimate SE df t.ratio p.value lower.CL upper.CL significance
(Medium-responders 1 Set) - (Medium-responders Baseline) -1.096 0.256 48.000 -4.279 0.003 -1.928 -0.264 **
(Medium-responders 1 Set) - (High-responders Baseline) -1.005 0.297 63.536 -3.379 0.032 -1.960 -0.050 **
(Medium-responders 1 Set) - (Low-responders 1 Set) -0.990 0.297 63.536 -3.329 0.036 -1.945 -0.035 **
(Low-responders 4 Sets) - (Medium-responders 1 Set) 0.986 0.297 63.536 3.315 0.038 0.031 1.941 **
(Medium-responders 4 Sets) - (Medium-responders Baseline) -0.932 0.256 48.000 -3.636 0.018 -1.763 -0.100 **
$LV38
LV38
contrast estimate SE df t.ratio p.value lower.CL upper.CL significance
(High-responders 4 Sets) - (High-responders Baseline) 0.403 0.095 48 4.247 0.003 0.095 0.711 **

Figures

LVs of interest

Figures for the AUC and FDR

  • Using this code only to adjust the names from the PLIER matrix
# Need to change names from the plierResults$U. 
rownames(plierResults$U)[rownames(plierResults$U) == "MetaMEx_AcuteResistance"] <- "MetaMex: Acute resistance exercise"
rownames(plierResults$U)[rownames(plierResults$U) == "REACTOME_METABOLISM_OF_PROTEINS"] <- "REACTOME: Metabolism of proteins"
rownames(plierResults$U)[rownames(plierResults$U) == "REACTOME_HIV_INFECTION"] <- "REACTOME: HIV infection"
rownames(plierResults$U)[rownames(plierResults$U) == "KEGG_PROTEASOME"] <- "KEGG: Proteasome"
rownames(plierResults$U)[rownames(plierResults$U) == "REACTOME_IMMUNE_SYSTEM"] <- "REACTOME: Immune system"
rownames(plierResults$U)[rownames(plierResults$U) == "MIPS_PA700_20S_PA28_COMPLEX"] <- "MIPS: PA700 20S PA28 complex"
rownames(plierResults$U)[rownames(plierResults$U) == "REACTOME_CLASS_I_MHC_MEDIATED_ANTIGEN_PROCESSING_PRESENTATION"] <- "REACTOME: Class I MHC mediated antigen processing presentation"
rownames(plierResults$U)[rownames(plierResults$U) == "MIPS_SPLICEOSOME"] <- "MIPS: Spliceosome"
rownames(plierResults$U)[rownames(plierResults$U) == "REACTOME_TRANSLATION"] <- "REACTOME: Translation"
rownames(plierResults$U)[rownames(plierResults$U) == "KEGG_SPLICEOSOME"] <- "KEGG: Spliceosome"
rownames(plierResults$U)[rownames(plierResults$U) == "REACTOME_METABOLISM_OF_LIPIDS_AND_LIPOPROTEINS"] <- "REACTOME: Metabolism of lipids and lipoproteins"
rownames(plierResults$U)[rownames(plierResults$U) == "REACTOME_METABOLISM_OF_AMINO_ACIDS_AND_DERIVATIVES"] <- "REACTOME: Metabolism of amino acids and derivatives"

Heatmap

# Need to source the custom PLIER function. KL and ML have worked on that.
source("Functions/PLIER_functions.R")

simpledf <- design_data_temp |> 
  dplyr::select(SubjectID, Sex, Group, TimePoint, Set) |> 
  arrange(Group, TimePoint, Set, Sex) |> 
  mutate(time = factor(paste(TimePoint, Set, sep = "_"))) |> 
  dplyr::select(Group, Sex, time, SubjectID) |> 
  rename("Groups"  = Group, "Sex" = Sex, "Condition" = time) |> 
  mutate(Sex = if_else(Sex == 1, "Female", "Male")) |> 
  mutate(Condition = factor(Condition, levels = c("pre_1", "post_1", "post_4"), labels = c("Baseline", "1 set", "4 sets"))) |> 
  mutate(Groups = factor(Groups, levels = c("L1_L4", "L1_H4", "H1_H4"), labels = c("Low-responders", "Medium-responders", "High-responders"))) |> 
  arrange(Groups, Condition) |> 
  column_to_rownames("SubjectID")

# Copy data.frame with normalized gene expression
newyN <- yN[, rownames(simpledf)]

# I prefer to personalize the annotation color here.
# The colors match all other graphs.
ann_colors = list(Condition = c(Baseline = "white", `1 set` = "#e2e2e2", `4 sets` = "#a8a8a8"),
    Sex = c(Female = "red", Male = "blue"),
    Groups = c(`Low-responders` = "#fdbeb8", `Medium-responders` = "#9bdaa3", `High-responders` = "#b3cefe"),
    present = c(inPathway = "black", notInPathway = "beige"))

# Run the custom pheatmap
LV_heatmaps <- LVs |> 
  map(\(x) gsub("LV", "", x)) |> 
  map(\(x) x = as.numeric(x)) |> 
  walk(\(x) ml1.customplotTopZ(plierRes = plierResults, data = newyN[cm.genes, ], priorMat = allPaths[cm.genes, ], top = 50, index = x, allLVs = T, annotation_col = simpledf, width = 30, height = 60, gaps_col = c(27,54), annotation_colors = ann_colors, annotation_legend = T, main = paste0("LV", x))) |> 
  set_names(LVs)